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Executive  Summary 


This  report  discusses  the  research  on  “Statistical  Signal  Analysis  Using  Wavelets” 
performed  by  the  Statistical  Sciences  Division  of  MathSoft,  Inc.  for  contract  N00014- 
93-C-0106  with  the  Office  of  Naval  Research.  The  overall  goals  of  the  research  are: 

•  Application  of  wavelets  and  related  transforms  to  data  analysis. 

•  Using  wavelets  as  the  basis  for  statistical  problems,  such  as  signal  extraction, 
spectral  density  estimation,  isotonic  regression,  and  classification. 

•  Research  into  new  transforms  and  algorithms  tailored  to  meet  the  needs  of 
data  analysis  and  statistical  estimation. 


The  first  year  of  reseaxch  focused  on  six  specific  axeas: 

1.  Applications  of  wavelets  to  data  of  interest  to  the  Navy; 

2.  Noise  removal  using  wavelets,  wavelet  packets,  and  cosine  packets; 

3.  Investigation  into  new  wavelet  transforms  which  are  outlier  resistant  and  edge 
preserving; 

4.  Development  of  a  framework  and  tools  for  the  “wavelet  approach”  towards 
analysis  of  signals,  images,  and  other  data; 

5.  Exploration  of  the  use  of  wavelets  as  a  dimension  reduction  tool  for  statistical 
analysis  of  very  large  data  sets; 

6.  Development  of  algorithms  for  wavelet  analysis. 


Future  research  will  extend  the  results  obtained  in  these  areas.  In  addition,  research 
in  the  following  areas  will  be  pursued,  time  permitting:  analysis  of  wavelet  coeffi¬ 
cients  for  fractional  Brownian  motion,  and  simulation,  bootstrapping,  and  modeling 
of  non-Gaussian  processes. 
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1  INTRODUCTION 


In  the  past  few  years,  wavel  ts  have  evoived  from  an  interesting  mathematical  dis¬ 
covery  to  a  valuable  technique  with  a  wide  variety  of  applications.  Wavelet  research 
has  synthesized  a  number  of  related  ideas  into  a  coherent  set  of  tools  and  method¬ 
ology  for  analysis  of  signals,  images  and  other  technical  data.  Research  is  leading 
to  new  wavelet-like  methods,  such  as  wavelet  packets  or  outlier-resistant  wavelets. 
The  development  of  software  toolkits  are  bringing  these  methods  to  the  hands  of 
scientists  and  engineers,  including  those  without  wavelet  expertise. 

Our  research  for  ONR  contract  N00014-93-C-0106  has  focused  on: 

•  Application  of  wavelets  and  related  transforms  to  data  analysis. 

•  Using  wavelets  as  the  basis  for  statistical  problems,  such  as  signal  extraction, 
spectral  density  estimation,  isotonic  regression,  and  classification. 

•  Research  into  new  transforms  and  algorithms  tailored  to  meet  the  needs  of 
data  analysis  and  statistical  estimation. 

Our  research  in  each  of  these  areas  is  guided  by  data  analysis  of  time  series  of  interest 
to  the  Navy. 

In  section  2,  we  summarize  the  main  results  from  the  first  year  of  research. 
Section  3  gives  a  list  of  all  papers  and  talks  produced  during  the  year.  Future 
research  plans  are  discussed  in  section  4.  We  also  include  an  appendix  with  plots 
giving  examples  and  illustrations. 

2  SUMMARY  OF  RESULTS 

We  are  pursuing  research  in  six  related  areas: 

1.  Applications  of  wavelets  to  data  of  interest  to  the  Navy; 

2.  Noise  removal  using  wavelets,  wavelet  packets,  and  cosine  packets; 

3.  Investigation  into  new  wavelet  transforms  which  are  outlier  resistant  ajid  edge 
preserving; 

4.  Development  of  a  framework  and  tools  for  the  “wavelet  approach”  towards 
analysis  of  signals,  images,  and  other  data; 

5.  Exploration  of  the  use  of  wavelets  as  a  dimension  reduction  tool  for  statistical 
analysis  of  very  large  data  sets; 

6.  Development  of  algorithms  for  wavelet  analysis. 

These  are  discussed  in  more  detail  below. 
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2.1  DATA  ANALYSIS  WITH  WAVELETS 

Applications  of  interest  to  the  Navy  include: 

•  Robust  wavelet  de-noising  of  radar  glint  noise. 

•  Time-frequency  analysis  of  underwater  acoustic  signals. 

•  Wavelet  compression  applied  to  a  number  of  signals  and  images. 

•  Fast  classification  of  transients  in  low  frequency  sinusoidal  data. 

Some  of  these  applications  are  discussed  in  more  detail  in  the  following  sections. 

We  are  currently  writing  up  a  series  of  reports  on  applications  of  wavelets  to 
data  analysis  problems  [BG95c,  BG95d,  BG95e].  Our  analysis  of  Navy  data  is  also 
closely  related  to  our  work  in  developing  the  “wavelet  approach”  to  data  analysis: 
see  section  2.4. 

2.2  NOISE  REMOVAL  WITH  WAVELETS 

In  the  past  three  years,  wavelet  de-noising  research  has  received  intense  activity  since 
the  initial  development  by  Donoho  and  Johstone.  This  research  includes  methods 
for  the  following  topics. 

Signal  recovery:  The  radar  glint  noise  example  is  an  example  of  signal  recovery  from 
noisy  data.  The  underlying  model  for  our  data  y,-  is  yi  —  fi  -f  e,-  where  /,•  is 
the  unknown  signal  or  image  and  Cj  is  noise.  Wavelets  provide  a  way  to  obtain 
an  estimate  /,■  while  making  a  minimum  of  assumptions  about  the  nature  of 
ft  and  Ci.  Other  applications  have  been  explored  by  [Bri94,  CM92,  DJ92a, 
DJ92b,  MH92,  Don93,  MZ93,  cA94,  HK94,  IK93,  PBBA94,  RBL94,  SS94, 
TH94,  Tew94]. 

Inverse  Problems  (Parameter  Estimation) :  Many  interesting  problems  having  to  do 
with  noisy  data  involve  indirect  measurements  yi  —  (/<’/)(<,)-{-€,•  where  we  want 
to  estimate  /  (so  called  inverse  problems).  Examples  of  the  transform  K  in¬ 
clude  the  Fourier  transformation  (magnetic  resonance  imaging),  the  Laplace 
transformation  (fluorescence  spectroscopy),  the  Radon  transformation  (tomog¬ 
raphy  problem)  and  various  deconvolution  problems  (gravity  anomalies,  in¬ 
frared  spectroscopy,  extragalactic  astronomy).  See  [Don93,  Wic94,  BFCLB94, 
MW94]  for  applications  of  wavelets  to  inverse  problems. 

Signal  Detection  and  Classification:  Closely  related  to  the  problem  of  signal  recov¬ 
ery  is  signal  detection  and  classification.  The  aim  is  to  identify  a  particular 
signal  from  background  noise  and  other  signals.  Wavelet  based  signal  detection 
has  been  explored  by  [FM92,  LKW92,  KDP92,  YRKS92,  Car94,  DMWJ94, 
LAJ94]. 
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Image  and  Image  Array  Reconstruction  and  Enhancement:  Wavelet  de-noising  method¬ 
ology  can  be  extended  to  two  and  higher  dimensional  problems  in  a  straight¬ 
forward  manner.  Some  applications  in  which  wavelets  have  been  used  to  clean 
noisy  images  are  given  in  [DJL92,  SFAH92,  MH92,  Don93,  WRM+94]. 

Density  estimation:  Wavelets  can  be  used  to  recover  probability  densities  [JKP92] 
and  spectral  densities  [Mou93,  Gao93b]  from  noisy  data  (see  section  2.2.5). 

In  year  two  research,  we  have  studied 

1.  Wavelet  de-noising  using  non-decimated  transforms. 

2.  Variance  estimation  for  wavelet  de-noising. 

3.  Selecting  the  threshold  using  cross-validation. 

4.  Applications  of  wavelets  to  isotonic  regression. 

5.  Spectral  density  estimation  using  wavelets. 

6.  Diagnostics  for  assessing  wavelet  de-noising. 

This  are  discussed  in  more  detail  below. 

2.2.1  Wavelet  De-Noising  using  Non-Decimated  Transforms 

The  non-decimated  discrete  wavelet  transform  is  a  non-orthogonal  variant  of  the 
classical  DWT.  With  the  non-decimated  DWT,  starting  with  n  sample  signal  val¬ 
ues  you  end  up  with  ( J  -|- 1)  x  n  coefficients.  Unlike  the  classical  DWT,  which  has 
fewer  coefficients  at  coarser  scales,  each  scale  for  the  non-decimated  DWT  has  n 
coefficients.  The  non-decimated  DWT  is  over-sampled  at  coarse  scales.  This  over- 
sampling  can  enhance  the  visual  displays  and  lead  to  advantages  in  certain  problems, 
including  better  spatial  resolution,  less  loss  of  information,  and  translation  invari¬ 
ance  with  respect  to  the  signal.  Mallat  and  co-workers  used  the  non-decimated 
DWT  for  detecting  singularities  and  de-noising  signals  [MH92,  MZ92]. 

Instead  of  using  an  orthogonal  discrete  wavelet  transform,  the  Donoho  and  John¬ 
stone  WaveShrink  algorithm  can  be  applied  to  the  non-decimated  wavelet  transform. 
Figure  1  illustrates  this  using  a  noisy  sinusoid  with  a  discontinuity.  The  WaveShrink 
estimate  based  on  the  orthogonal  haar  wavelet  transform  is  very  blocky,  correspond¬ 
ing  to  the  discontinuous  nature  of  the  haar  wavelet.  By  contrast,  the  WaveShrink 
estimate  bcised  on  based  on  the  non-decimated  haar  wavelet  estimate  is  smooth, 
due  to  the  “spatial  averaging”  in  the  reconstruction.  There  is,  however,  a  cost  to 
smoothness:  the  non-decimated  estimate  blurs  the  discontinuities  at  the  two  jumps. 
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Preliminary  results  indicate  that  de-noising  with  non-decimated  wavelets  offers 
significant  improvements  over  the  original  WaveShrink  algorithm.  We  plan  to  con¬ 
tinue  our  research  into  de-noising  with  non-decimated  wavelets,  through  both  an 
empirical  study  and  the  development  of  software  for  use  by  Navy  researchers  at  the 
China  Lake  facility. 

2.2.2  Variance  Estimation  for  Wavelet  De-Noising 

Based  on  ideas  proposed  by  David  Brillinger  [Bri94],  we  have  developed  a  gen¬ 
eral  method  for  estimating  the  variance  of  WaveShrink  estimate.  Let  ft  be  the 
WaveShrink  estimate.  Since  the  wavelet  transform  is  linear  and  the  wavelet  coeffi¬ 
cients  Wj^k  are  approximately  uncorrelated,  asymptotic  theory  implies  that 

Var(/t)  Ri  Var(rwj,A:)  (1) 

i,k 

where 

Wj,k  =  Sx{wj,k)wj,kUj,k{t)\\  (2) 

In  figure  2,  we  apply  this  formula  to  compute  a  confidence  interval  for  the  flow 
of  the  Nile  River.  The  top  plot  is  the  annual  Nile  River  flow  at  Aswan  from  1875  to 
1970,  the  middle  plot  is  the  WaveShrink  using  Haar  wavelet  with  an  approximated 
95%  confidence  band  and  the  bottom  plot  is  the  WaveShrink  using  bsl.3  wavelet 
with  an  approximated  95%  confidence  band.  The  confidence  interval  clearly  shows 
a  significant  drop  around  the  time  of  the  building  of  the  Aswan  dam  in  1899-1902. 

To  make  (1)  more  generally  useful,  we  plan  to  investigate  several  extensions. 
These  include  the  development  of  a  fast  algorithm  to  compute  \\ipj,k{t)\\  and  use  of 
the  bootstrap  to  get  a  more  realistic  estimate  of  Var(thj,fc). 

2.2.3  Cross  Validation  Selection  of  the  Threshold 

Initial  simulations  have  shown  that  the  cross  validation  technique  is  promising  in 
selecting  thresholds  in  the  WaveShrink  algorithm.  The  basic  idea  is: 

[1]  for  each  set  of  thresholds,  apply  WaveShrink  to  part  of  the  data  and  compute 
the  mean-square  error  between  the  WaveShrink  estimate  and  the  other  part 
of  the  data; 

[2]  find  the  set  of  thresholds  which  minimizes  the  mean-square  errors. 

In  figure  3,  we  apply  cross  validation  to  select  the  threshold  for  WaveShrink 
applied  to  the  synthetic  “doppler”  signal  corrupted  by  Gaussian  noise.  The  cross- 
validation  threshold  is  compared  with  the  “universal”  threshold,  “minimax”  thresh¬ 
old  and  the  “optimal”  threshold  (the  threshold  minimizes  the  mean-square  error 
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(MSE)  between  WaveShrink  estimate  and  the  true  signal,  this  is  only  possible  for 
simulation).  Cross  validation  consistently  leads  to  lower  MSE  than  the  “universal” 
and  “minimax”  thresholds  for  this  example. 

In  future  research,  we  will  to  explore  using  non-linear  optimization  techniques 
to  estimate  multiple  thresholds  with  cross  validation.  We  will  also  explore  other 
cross  validation  methods,  such  as  using  decimation  by  powers  of  2-’  instead  of  just 
decimation  by  2  (initial  investigations  have  not  shown  this  to  be  crucial,  but  we 
suspect  this  may  not  be  the  case  for  “real  world”  problems).  We  will  also  develop  a 
theoretical  justification  for  cross-validation  in  the  WaveSrhink  context. 

2.2.4  Wavelets  and  Isotonic  Regression 
Consider  the  isotonic  regression  model: 

Vi  =  f{U)  +  Zi 

where  /  is  a  decreasing  function  and  {zi}  are  assumed  to  be  a  stationary  Gaussian 
process  with  mean  zero  and  variance  cr^.  We  propose  a  simple  thresholding  procedure 
beised  on  the  fact  that  the  wavelet  coefficients  for  /,  under  Haar  wavelet,  are  non¬ 
negative.  We  show  that  our  estimator  is  competitive  with  the  Grenander  estimator 
both  theoretically  and  numerically  (in  the  sense  of  mean-square-error).  Figure  4 
displays  a  synthetic  decreasing  curve  (top  left),  the  same  curve  plus  Gaussian  white 
noise  (top  right),  the  Grenander  estimate  (bottom  left)  and  the  wavelet  estimate 
(bottom  right).  The  wavelet  estimate  has  lower  mean  square  error  and  does  a  better 
job  of  preserving  the  jumps. 

Details  of  the  wavelet  estimation  procedure  are  given  in  [Gao95b]:  see  section  3. 

A  limitation  of  the  current  procedure  is  the  restriction  to  the  use  of  the  Haar 
wavelet,  which  is  the  only  wavelet  which  preserves  the  monotonicity.  We  will  inves¬ 
tigate  the  existence  of  a  new  class  of  smooth  wavelets  which  preserve  monotonicity. 

2.2.5  Wavelet  Based  Spectral  Density  Estimation 

A  technique  for  spectral  density  estimation  based  on  wavelet  decomposition  of  the 
periodogram  and  reconstruction  of  the  spectrum  was  developed  by  Hong- Ye  Gao 
[Gao93b,  Gao93a].  This  technique  is  especially  valuable  for  processes  with  non¬ 
smooth  features  in  the  spectrum,  such  as  sharp  peaks. 

A  simulated  autoregressive  series  is  displayed  in  figure  5.  This  series  has  sharp 
peaks  in  its  spectral  density  function.  Figure  6  compares  wavelet  shrinkage  estima¬ 
tion  applied  to  this  series  to  a  more  traditional  non-parametric  spectral  estimator 
based  on  a  kernel  smoother.  The  traditional  nonparametric  estimator  based  on  a 
triangular  spectral  smoothing  window  tends  to  oversmooth  the  peaks.  Using  shorter 
span  smoothing  windows  would  better  preserve  the  peaks,  but  would  result  an  un¬ 
necessarily  rough  estimate  elsewhere.  The  WaveShrink  estimator  is  equivalent  to 
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using  a  variable  bandwidth  smoother.  It  preserves  the  peaks  while  producing  a 
smooth  estimate  elsewhere. 

A  paper  describing  Dr.  Gao’s  work  has  been  accepted  by  the  Journal  of  Time 
Series  Analysis  pending  revisions  [Gao95a]  (see  section  3).  These  revisions  are  being 
done  through  the  support  of  this  contract. 

2.2.6  Diagnostics  for  Assessing  Wavelet  De-Noising 

The  theory  for  WaveShrink  is  based  on  the  Gaussian  white  noise  model.  When  this 
model  is  inappropriate  -  e.g.,  when  the  data  has  outliers  -  then  the  “traditional” 
WaveShrink  estimators  may  be  inappropriate. 

To  help  assess  the  WaveShrink  fit,  we  have  explored  the  use  of  simple  diagnostics 
plots.  Figure  7  displays  a  synthetic  bumps  signal,  the  bumps  signal  corrupted  by 
Gaussian  noise,  and  the  WaveShrink  estimate  of  the  bumps  signal.  How  should  we 
assess  the  WaveShrink  estimate?  Figure  8  gives  four  “views”  of  the  WaveShrink 
estimate  for  the  bumps  signal:  (1)  The  decomposition  of  the  data  into  signal  plus 
noise,  (2)  Boxplots  of  the  DWT  coefficients  for  the  original  data  with  the  WaveShrink 
thresholds  superimposed  as  horizontal  lines,  (3)  The  DWT  of  the  signal,  and  (4) 
A  barplot  showing  the  decomposition  of  the  energy  of  the  data  into  the  energy 
attributable  to  signal  and  residual  energy. 

We  can  also  use  plots  to  examine  the  residuals  from  the  figure.  Figure  9  dis¬ 
plays  a  series  of  views  for  the  residuals  from  the  WaveShrink  estimate  of  the  bumps 
signal.  The  diagnostic  plots  indicate  that  the  peaks  are  oversmoothed.  As  a  re¬ 
sult,  the  distribution  of  the  residual  component  is  skewed  toward  high  values.  The 
oversmoothing  also  leads  to  significant  autocorrelation  in  the  residuals. 

We  will  continue  to  explore  the  use  of  visual  diagnostics  to  help  guide  the  user 
when  to  use  WaveShrink,  and  how  to  improve  the  WaveShrink  fit. 

2.3  ROBUST  AND  NONLINEAR  WAVELETS 

The  aim  of  this  research  is  to  investigate  wavelet  methods  which  are  robust  towards 
outliers.  This  line  of  research  was  motivated  by  problems  encountered  in  application 
of  wavelets  to  Navy  data  sets.  It  promises  to  be  of  considerable  practical  importance, 
and  represents  an  exciting  area  of  innovative  research. 

'Some  outlier  resistant  and  edge  preserving  wavelets  we  developed  include: 

1.  Wavelets  with  robust  smoother /cleaner.  The  usual  wavelet  transform  is  com¬ 
bined  with  a  robust  smoother/ cleaner  to  remove  outliers. 

2.  Wavelet  based  minimum  entropy  segmentation.  The  wavelet  transform  is  used 
to  segment  a  noisy  signal  by  successively  identifying  and  removing  edges  and 
other  discontinuities. 
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These  are  discussed  in  more  detail  below. 

Research  in  this  area  is  still  continuing:  see  sections  3.6  and  4. 

2.3.1  Wavelets  with  Robust  Smoother/Cleaner 

The  presence  of  outliers  in  data  causes  problems  in  traditional  time  series  analysis 
techniques.  Outliers  can  seriously  distort  the  autocorrelation  function,  partial  auto¬ 
correlation  function,  spectral  density  function,  model  identification,  and  parameter 
estimates  for  models.  Outliers  can  also  cause  problems  with  methods  based  on  the 
wavelet  decomposition.  Wavelets  are  a  linear  transformation  of  the  data,  and  hence, 
outliers  have  unbounded  influence  on  the  wavelet  coefficients. 

The  goal  of  robust  smoother/ cleaner  wavelets  is  to  produce  a  fast  wavelet  decom¬ 
position  which  is  robust  towards  outliers.  Smoother-cleaner  wavelets  behave  like  the 
classical  L2  wavelet  transform  for  Gaussian  signals,  but  prevent  outliers  and  outlier 
patches  from  leaking  into  the  wavelet  coefficients  at  coarse  levels  (like  Li  wavelets). 
However,  in  contrast  to  the  wavelets,  algorithm  is  very  fast  with  computational 
complexity  0{n). 

The  basic  idea  of  robust  smoother/cleaner  wavelets  is  simple:  the  smooth  coef¬ 
ficients  are  preprocessed  with  a  fast  and  robust  smoother/cleaner.  The  procedure 
is  illustrated  in  figure  10.  As  usual,  we  start  with  a  set  of  wavelet  coefficients  5(0). 
Then,  for  each  multiresolution  level,  the  signal  is  decomposed  into  three  components: 

1.  A  set  of  robust  residuals  R(i  —  1),  given  by 

R{e  -l)  =  Sx  [S{£  -  1)  -  S{£  -  1)) 

where  6\  is  a  thresholder  function  and  5(^  —  1)  is  a  robust  smooth  of  5(£  —  1) 
(e.g.,  running  medians  of  5).  The  threshold  A  is  chosen  so  that  most  of  the 
robust  residuals  are  zero. 

2.  The  smooth  wavelet  coefficients  5(^)  obtained  by  applying  the  usual  low- 
pass/decimation  wavelet  filter  H  to  the  cleaned  smooth  coefficients  U{£—1)  = 

S{£-1)-Ri£-1). 

3.  The  detail  wavelet  coefficients  D{£)  obtained  by  applying  the  usual  high- 

'  pass/decimation  wavelet  filter  G  to  U{£  —  1). 

The  smoother/cleaner  wavelet  decomposition  can  be  used  for  de-noising  signals. 
As  an  example,  we  apply  it  to  the  radar  glint  noise  signal.  The  original  noisy  signal, 
which  is  the  angle  of  the  target  in  degrees,  is  displayed  in  Figure  11(a)  The  true 
signal  is  a  low-frequency  oscillation  about  0°.  The  signal  contains  a  number  of  glint 
spikes,  causing  the  apparent  signal  to  be  different  from  the  true  signal  by  as  much 
as  150°. 
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Figure  11(b)  compares  denoising  with  linear  shrinkage  of  wavelets  (dashed  line) 
to  denoising  with  WaveShrink  combined  with  robust  smoother-cleaner  wavelets 
(solid  line).  The  lineax  shrinkage  is  based  on  annihilating  all  detail  coefficients  of  the 
classical  wavelet  transform  at  levels  I  —  1,2, 3,4.  While  linear  shrinkage  estimate 
is  smooth,  it  is  still  somewhat  sensitive  to  the  glint  spikes.  By  contrast,  the  clean 
and  repeat  procedure  is  quite  resistant  to  the  glint  spikes  but  effectively  tracks  the 
sudden  changes  in  the  series. 

This  work  led  to  publication  of  two  papers  [BDGM94a,  BDGM94b]:  see  section  3. 
Refer  to  these  papers  for  details  about  the  smoother-cleaner  wavelets. 

2.3.2  Segmented  Wavelet  Bases 

Another  way  to  generalize  wavelet  bases  is  through  segmentation.  Segmented  wavelet 
bases  are  obtained  by  piecing  together  different  wavelet  bases  over  time  and  space. 
Segmented  bases  are  particularly  well  suited  for  adaptation  to  discontinuities  (edges, 
cusps,  etc.)  and  more  general  change  points  (e.g.,  changes  in  the  stochastic  prop¬ 
erties  of  the  signal).  The  segmentation  is  done  by  optimizing  some  criterion  (e.g., 
entropy)  in  a  similar  manner  cis  in  the  search  for  a  “best  basis”. 

In  work  supported  by  this  contract,  the  paper  “On  Minimum  Entropy  Segmen¬ 
tation”  will  appear  in  a  future  volume  of  Wavelet  Analysis  and  Its  Applications, 
edited  by  Charles  Chui  [Don95]  (see  section  3). 

2.4  A  FRAMEWORK  FOR  THE  “WAVELET  APPROACH” 

In  the  late  1960’s  and  early  1970’s,  John  Tukey  developed  and  popularized  a  phi¬ 
losophy  for  analyzing  data  called  “exploratory  data  analysis”  [Tuk75].  Since  that 
time,  this  philosophy  has  expanded  and  matured  into  a  complete  collection  of  tools 
and  techniques  for  understanding  and  extracting  information  from  data.  Our  belief 
is  that  wavelet  analysis  will  follow  a  similar  evolution.  What  started  as  an  interest¬ 
ing  mathematical  breakthrough  has  evolved  into  a  coherent  new  way  of  analyzing 
signals,  images,  and  other  data. 

Many  of  the  ideas  behind  wavelets  are  drawn  from  other  domains,  such  as  sub¬ 
band  filtering,  approximation  theory,  signal  processing,  and  image  processing.  What 
is  significant  about  wavelets  is  the  development  of  a  coherent  framework,  which  pro¬ 
vides  a  new  paxadigm  for  analyzing  data.  We  are  now  at  the  stage  where  it  is 
possible  to  formulate  this  framework  in  an  organized  manner  which  is  accessible  to 
the  broader  community  of  scientists,  engineers,  and  researchers,  and  not  just  wavelet 
experts. 

One  aspect  of  our  research  is  to  develop  this  framework,  fueling  the  transition  of 
wavelets  from  a  specialized  technique  to  a  broadly  used  methodology.  Our  efforts  in 
this  direction  have  resulted  in  the  writing  of  an  article  for  publication  in  the  IEEE 
Spectrum  [BDG95]:  see  section  3. 
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Another  aspect  of  research  related  to  developing  the  “wavelets  approach”  is  the 
design  of  S+ WAVELETS  toolkit.  The  design  of  the  toolkit  is  closely  linked  to  our 
efforts  to  develop  a  coherent  framework  for  wavelet  analysis.  The  design  of  the 
toolkit  is  discussed  in  the  technical  report  [BG95b].  The  toolkit  design  and  the 
corresponding  technical  report  were  partially  supported  by  this  contract. 

We  will  continue  our  research  in  this  area,  both  in  regards  to  writing  new  papers 
and  developing  software  promoting  the  use  and  application  of  wavelets. 

2.4.1  S+WAVELETS  Toolkit 

S+WavelETS  is  an  object-oriented  toolkit  for  wavelet  analysis  of  time  series  and 
images  [BG94].  It  is  the  only  commercial  toolkit  for  wavelet  analysis  incorporated 
into  a  high-level  interactive  language  such  as  S-PLUS.  A  high-level  overview  of  the 
toolkit  is  given  in  figures  14  and  15. 

The  primary  objective  of  the  S-f  WAVELETS  toolkit  is  to  provide  a  mature  com¬ 
mercial  wavelet  analysis  product  to  the  engineering  and  scientific  research  com¬ 
munity.  The  commercial  software  for  S-f- WAVELETS  is  being  developed  under  the 
support  of  NASA  Phase  II  SBIR  Contract  No.  NAS13-587.  The  S-(- Wavelets 
toolkit  provides  important  support  for  our  research  program  with  ONR: 

•  It  provides  an  environment  in  which  we  can  rapidly  prototype  and  test  new 
wavelet  based  methods  which  we  develop  in  our  ONR  research.  This  includes 
not  only  testing  methods  to  determine  performance  on  axtificial  or  real  test 
data  sets,  but  also  evaluating  the  practical  implementations  on  real  data  sets 
of  interest  to  the  Navy. 

•  The  existence  of  our  NASA  SBIR  funded  S-f  WAVELETS  commercial  software 
development  work  means  that  we  can  provide  China  Lake  (Gary  Hewer  and 
colleagues)  with  specialized  software  for  wavelet  analysis. 

We  will  continue  to  develop  software  based  on  the  toolkit  for  use  in  our  research 
and  by  Navy  researchers  at  China  Lake.  In  particular,  we  will  develop  software 
for  efficient  wavelet  processing  of  signal  and  image  arrays  and  non-decimated  2-D 
wavelet  transforms. 

2.5  DIMENSION  REDUCTION 

Since  wavelets  are  very  effective  at  compacting  energy,  wavelets  are  good  tools  for 
dimension  reduction  in  problems  involving  very  large  data  sets.  The  dimension 
reduction  property  of  wavelets  is  important  for  applications  such  as  image  compres¬ 
sion,  factor  analysis,  and  numerical  analysis:  see  chapter  11  of  [Wic94]. 

We  are  investigating  the  use  of  wavelets  as  a  dimension  reduction  tool  for  the 
problem  of  classification  of  transient  acoustic  signals  (see  below).  In  future  research. 
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we  may  also  look  at  use  of  dimension  reduction  as  a  means  to  apply  other  statistical 
techniques  to  very  large  data  sets. 

2.5.1  Algorithms  for  Classification 

Following  the  work  by  Saito  and  Coifman  [SC94],  we  illustrate  the  basic  ideas  of 
using  wavelets  to  reduce  the  dimension  of  the  problem  of  classifying  transient  acous¬ 
tic  signals.  We  excerpt  three  classes  of  signals:  whale  clicks,  snapping  shrimp  and 
background  noise.  Figure  12  displays  three  samples  of  length  2048  from  the  three 
classes.  The  wavelet  clcissification  scheme  is  based  on  reducing  the  signals  of  di¬ 
mension  2048  to  vectors  of  length  6,  and  building  a  classification  tree  based  on  the 
length  6  vectors.  Figure  13  shows  the  resulting  classification  tree. 

Specifically,  this  tree  was  built  as  follows: 

[1]  25  samples  from  each  class  are  drawn  to  form  a  training  data  set. 

[2]  For  each  sample,  a  wavelet  packet  table  is  computed  and  then  an  energy  map 
is  obtained  by  squaring  the  wavelet  packet  table  and  rescaling  by  the  total 
energy  of  the  signal. 

[3]  The  energy  maps  from  each  class  are  then  combined  to  form  a  single  energy 
map  for  the  class. 

[4]  Distances  between  the  energy  maps  (totally  3  in  this  example)  are  computed 
and  combined  to  a  single  distance  measure.  An  optimal  basis  which  maximizes 
the  distance  is  derived.  The  optimal  basis  in  this  example  consists  of  crystals: 
w4.0,  w4.1,  w4.2,  w4.3,  w2.1,  wl.l.  L2  distance  is  used  in  this  example, 
other  discrimination  measures  include  relative  entropy  and  general  Lp  distance^ 
[SC94]. 

[5]  The  training  samples  are  transformed  to  the  optimal  basis  and  total  energy 
for  each  crystal  is  accumulated. 

[6]  A  classification  tree  is  grown  based  on  the  accumulated  energy  vectors  (of 
length  6  in  this  example)  and  is  displayed  in  Figure  13. 

25Q  samples  are  drawn  from  each  class  and  classified  by  the  classification  tree  (apply 
step  [5]  above).  The  misclassification  rate  is  13.33%. 

In  future  research,  we  will  investigate  some  of  the  following  improvements 

[1]  Shrink  the  energy  maps  to  eliminate  ambient  noise  and  enhance  features. 

[2]  Time-invariant  discriminant  measures,  such  as  rank  related  (e.g.  Kruskal- 
Wallis  [BD77])  statistics  for  comparing  the  component- wise  distributions  of 
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the  transform  coefficients,  should  be  used  for  time-invariant  signals  (like  whale 
clicks).  Since  these  discrimination  meaisures  may  not  be  additive,  new  optimal 
basis  selection  algorithms  need  to  be  developed. 

[3]  For  a  given  optimal  transform,  an  equally  challenge  issue  is  how  to  compress 
the  transformed  coefficients.  For  time-variant  signals,  taking  certain  number 
of  top  coefficients,  cis  suggested  by  Saito  and  Coifman,  is  a  straight  forward 
solution.  However,  when  signals  are  time- in  variant,  this  is  not  be  the  case: 
a  transient  may  occur  at  any  time  and  therefore  the  coefficients  associated 
with  the  transient  may  occur  at  any  time.  In  above  example,  we  have  used 
the  accumulated  energy.  There  are  other  options  we  can  explore,  such  as  top 
frequencies,  autocorrelations  and  moments. 

2.6  ALGORITHMS  FOR  WAVELET  ANALYSIS 

To  support  our  research  for  analysis  of  Navy  data,  we  have  spent  some  effort  in 
development  of  fundamental  algorithms  for  wavelet  analysis.  This  work  has  involved 
two  separate  projects: 

•  Investigation  in  the  fundamental  properties  of  certain  biorthogonal  wavelet 
functions. 

•  Development  of  algorithms  to  handle  signals  and  images  of  arbitrary  dimen¬ 
sions  with  a  vaxiety  of  boundary  conditions. 

This  work  is  described  in  more  detail  below. 

2.6.1  Fundamental  Properties  of  Biorthogonal  Wavelets 

In  development  of  the  S-|- WAVELETS  toolkit,  we  discovered  a  problem  in  graphing 
certain  biorthogonal  wavelets;  namely,  the  wavelets  in  Tables  8.2  and  8.3  of  the  book 
by  Daubechies  [Dau92]  with  (N,  N)  =  (2, 2),  (3, 1),  (3, 3),  and  (5, 5).  In  the  toolkit, 
these  are  referred  to  as  wavelets  bs2.2,  bs3.1,  bs3.3,  and  vs3.  It  turns  out  that 
these  wavelets  are  infinite  at  all  dyadics!  This  fact  may  have  important  implications 
for  analysis  using  these  wavelets,  since  it  implies  a  certain  underlying  instability. 

We  were  naturally  led  to  the  question:  if  the  wavelets  are  infinite  at  all  dyadics, 
how  do  we  graph  the  wavelets  to  get  an  understanding  of  their  shape?  The  usual 
constructions  for  evaluating  wavelet  functions,  described  by  Strang  [Str89],  do  not 
apply. 

In  work  partially  supported  by  this  contract,  we  developed  an  approach  for 
graphing  such  wavelets.  This  approach  formulates  an  eigenvalue  equation  based 
on  the  dilation  equation  evaluated  at  non-dyadic  points.  This  work  is  described 
in  a  paper  published  in  “Recent  Advances  in  Approximation  Theory,  Wavelets  and 
Applications”  [RBG94]. 
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2.6.2  Algorithms  for  Finite  Signals 

A  very  important  problem  in  wavelet  analysis  of  signals  and  images  is  the  ability  to 
handle  arbitrary  signal/image  samples  sizes  with  a  variety  of  boundary  conditions. 
In  on-going  research,  we  are  developing  a  suite  of  algorithms  to  address  this  problem. 

Most  published  research  in  wavelet  analysis  is  based  on  signals  and  images  which 
axe  of  length/ dimension  2"^.  Unfortunately,  particularly  for  images,  we  do  not  always 
have  the  luxury  to  restrict  our  sample  size  in  that  manner.  Also,  in  wavelet  analysis, 
the  treatment  of  the  boundaries  is  very  important.  For  computational  simplicity, 
it  is  often  assumed  the  image  or  signal  is  periodic.  This  can  lead,  however,  to 
very  serious  artifacts  which  cause  problems  in  both  statistical  and  data  compression 
applications. 

In  a  technical  report,  partially  supported  by  this  contract,  we  describe  a  suite 
of  algorithms  for  addressing  these  issues  [BG95a].  In  work  in  progress,  we  are 
extending  this  work,  further  improving  these  algorithms  and  developing  a  conceptual 
framework  for  understanding  finite  wavelet  operators  [BGR95]:  see  section  3. 

3  PAPERS  AND  TALKS 

3.1  Published  Papers 

1.  Smoothing  and  Robust  Wavelet  Analysis.  Andrew  G.  Bruce,  David  L. 
Donoho,  Hong- Ye  Gao,  and  R.  Douglas  Martin.  In  Proceedings  in  Computa¬ 
tional  Statistics  (invited  paper),  Vienna,  Austria,  August,  1994. 

In  a  series  of  papers,  Donoho  and  Johnstone  develop  a  powerful  theory  based 
on  wavelets  for  extracting  non- smooth  signals  from  noisy  data.  Several  non¬ 
linear  smoothing  algorithms  are  presented  which  provide  high  performance 
for  removing  Gaussian  noise  from  a  wide  range  of  spatially  inhomogeneous 
signals.  However,  like  other  methods  based  on  the  linear  wavelet  transform, 
these  algorithms  are  very  sensitive  to  certain  types  of  non-Gaussian  noise,  such 
as  outliers.  In  this  paper,  we  develop  outlier  resistant  wavelet  transforms. 
In  these  transforms,  outliers  and  outlier  patches  are  localized  to  just  a  few 
scales.  By  using  the  outlier  resistant  wavelet  transforms,  we  improve  upon  the 
Donoho  and  Johnstone  nonlinear  signal  extraction  methods.  The  outlier  resis¬ 
tant  wavelet  algorithms  are  included  with  the  S-f- WAVELETS  object-oriented 
toolkit  for  wavelet  analysis. 

2.  Nonlinear  and  Robust  Wavelet  Analysis.  Andrew  Bruce,  David  L. 
Donoho,  Hong- Ye  Gao,  and  R.  Douglas  Martin.  In  SPIE  Proceedings,  Wavelet 
Applications,  Orlando,  Florida,  April,  1994. 

See  above  for  abstract, 
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3.  Non-smooth  Wavelets:  Graphing  Functions  Unbounded  on  Every 
Interval.  David  Ragozin,  Andrew  Bruce,  and  Hong- Ye  Gao.  In  Recent  Ad¬ 
vances  in  Approximation  Theory,  Wavelets  and  Applications,  Kluwer,  1994. 

Several  wavelets  from  well  known  biorthogonal  families  are  shown  to  be  un¬ 
bounded  on  every  interval.  One,  in  fact,  is  shown  to  be  infinite  at  each  dyadic 
rational.  Not  withstanding  these  facts,  we  show  how  to  compute  exact  values 
for  these  wavelets  at  many  points  and  thus  achieve  exact  pictures  for  these 
functions. 


3.2  Papers  Accepted  for  Publication 

1.  Ideal  Denoising  in  an  Orthonormal  Basis  Selected  from  a  Library 
of  Bases  David  L.  Donoho  and  Iain  M.  Johnstone.  To  appear  in  Comptes 
Rendus  de  V Academic  de  Science  Paris. 

Suppose  we  have  observations  y,  =  s,  Zi,  i  =  1,  . . .,  n,  where  (s,)  is  signal 
and  {zi)  is  i.i.d.  Gaussian  white  noise.  Suppose  we  have  avilable  a  library  C  of 
orthogonal  bases,  such  as  the  Wavelet  Packet  ba.ses  or  the  Cosine  Packet  bases 
of  Coifman  and  Meyer.  We  wish  to  select,  adpatively  based  on  the  noisy  data 
(?/,),  a  basis  in  which  best  to  recover  the  signal  (“de- noising”).  Let  Mn  be  the 
total  number  of  distinct  vectors  occurring  among  all  bcises  in  the  library  and 
let  tn  =  y^21og(M„).  (For  wavelet  packets,  Mn  =  nlog2(n).) 

Let  y[B\  denote  the  original  data  y  transformed  into  the  Basis  B.  Choose 
A  >  8  and  set  Lambdon  =  (A  •  (1  -f-  tn)Y.  Define  the  entropy  functional 

=  Z!niin(r/f[H],A2). 
i 

Let  B  be  the  best  orthogonal  basis  according  to  this  entropy: 

B  =  aTgimnSx{y,B). 

Define  the  hard-threshold  nonlinearity  Tit{y)  =  In  the  empirical  best 

basis,  apply  hard-thresholding  with  threshold  t  = 

'  s*[B]  =  v^{yi[B]). 

Theorem:  With  probability  exceeding  TTn  =  1  —  e/Mn, 

ll-s*  -  •s|l2  <  (1  -  8/^)""^  •  An  •  mmF||Js  -  s||^. 

Here  the  minimum  is  over  all  ideal  procedures  working  in  all  vases  of  the 
library,  i.e.  in  basis  B,  sg  is  just  i/i[H]l{|s;[e]>i}. 
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In  short,  the  basis-adaptive  estimator  achieves  a  loss  within  a  logarithmic 
factor  of  all  the  ideal  risk  which  would  be  achievable  if  one  had  available  an 
oracle  which  would  supply  perfect  information  about  the  ideal  basis  in  which 
to  de-noise,  and  also  about  which  coordinates  were  large  or  small. 

The  result  extends  in  obvious  ways  to  more  general  orthogonal  basis  libraries, 
basically  to  any  libraries  constructed  from  an  at-most  polynomially-growing 
number  of  coefficient  functionals.  Parallel  results  can  be  developed  for  closely 
related  entropies. 

2.  On  Minimum  Entropy  Segmentation.  David  L.  Donoho.  To  appear  in 
Wavelet  Analysis  and  Its  Applications,  edited  by  Charles  Chui. 

A  segmented  multiresolution  analyses  of  [0, 1]  is  described.  Such  multiresolu¬ 
tion  analyses  lead  to  segmented  wavelet  bases  which  are  adapted  to  disconti¬ 
nuities,  cusps,  etc.,  at  a  given  location  r  G  [0,1].  The  approach  emphasizes 
the  idea  of  average-interpolation  -  synthesizing  a  smooth  function  on  the  line 
having  prescribed  boxcar  averages.  This  particular  approach  leads  to  methods 
with  subpixel  resolution  and  to  wavelet  transforms  with  the  advantage  that, 
for  a  signal  of  length  n,  all  n  pixel- level  segmented  wavelet  transforms  can  be 
computed  simultaneously  in  a  total  time  and  space  which  axe  both  0{n  log(n)). 

The  search  for  a  segmented  wavelet  basis  is  considered  which,  among  all  such 
segmented  bases,  minimizes  the  “entropy”  of  the  resulting  coefficients.  Fast 
access  to  all  segmentations  enables  fast  search  for  a  best  segmentation. 

When  the  “entropy”  is  Stein’s  Unbiased  Risk  Estimate,  one  obtains  a  new 
method  of  edge-preserving  de-noising.  When  the  “entropy”  is  the  £^-energy, 
one  obtains  a  new  multi-resolution  edge  detector,  which  works  not  only  for 
step  discontinuities  but  also  for  cusp  and  higher-order  discontinuities,  and  in 
a  near-optimal  fashion  in  the  presence  of  noise. 

An  iterative  approach  is  also  described,  Segmentation  Pursuit,  for  identifying 
edges  by  the  fast  segmentation  algorithm  and  removing  them  from  the  data. 

3.3  Papers  Accepted  Pending  Review 

1.  Wavelet  Analysis  Tools  Become  Widely  Available.  Andrew  G.  Bruce, 
David  L.  Donoho,  and  Hong- Ye  Gao.  Under  review  at  the  IEEE  Spectrum. 

This  is  an  introductory  article  intended  to  promote  the  use  and  application  of 
wavelets  and  wavelet  software.  The  article  gives  a  non-mathematical  review 
of  wavelet  analysis,  complete  with  color  graphics  and  informative  sidebars. 
Several  applications  of  wavelets  are  given,  including  data  compression,  noise 
removal,  and  development  of  fast  algorithms.  The  discussion  includes  the 
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latest  developments  in  wavelets,  including  time-frequency  decompositions  such 
as  wavelet  packets  and  cosine  packets. 

The  main  focus  of  the  article  is  to  give  a  description  of  the  “wavelets  toolkit”, 
which  is  the  collection  of  software  tools  needed  to  perform  a  wavelet  analysis. 
A  wide  range  of  commercial  and  public  domain  software  is  reviewed.  The 
advantages  and  disadvantages  of  different  wavelet  computing  environments 
are  discussed. 

2.  Choice  of  Thresholds  for  Wavelet  Shrinkage  Estimate  of  the  Spec¬ 
trum.  Hong- Ye  Gao.  Under  review  at  the  Journal  of  Time  Series  Analysis. 

We  study  the  problem  of  estimating  the  log  spectrum  of  a  stationary  Gaussian 
time  series  by  thresholding  the  empirical  wavelet  coefficients.  We  propose  the 
use  of  thresholds  depending  on  sample  size  n,  wavelet  basis  and  resolution 
level  j.  At  fine  resolution  levels  {j  =  1,2,...),  we  propose 

~  OCj  log  Tl , 

where  {oj}  are  level- dependent  constants  and  at  coarse  levels  {j  >•  1) 

The  purpose  of  this  thresholding  level  is  to  make  the  reconstructed  log-spectrum 
as  nearly  noise-free  as  possible.  In  addition  to  being  pleasant  from  a  visual 
point  of  view,  the  noise-free  character  leads  to  attractive  theoretical  properties 
over  a  wide  range  of  smoothness  assumptions.  Previous  proposals  set  much 
smaller  thresholds  and  did  not  enjoy  these  properties. 

3.4  Technical  Reports 

1.  Wavelet  and  Isotonic  Regression.  Hong- Ye  Gao. 

Consider  the  model: 

Vi  =  f{ti)  +  Zi 

where  /  is  a  decreasing  function  and  {zi}  are  assumed  to  be  a  stationary 
Gaussian  process  with  mean  zero  and  variance  cr^.  We  propose  a  simple  thresh¬ 
olding  procedure  based  on  the  fact  that  the  wavelet  coefficients  for  /,  under 
Haar  basis,  are  non-negative.  We  show  that  our  estimator  is  competitive  with 
the  Grenander  estimator  both  theoretically  and  numerically  (in  the  sense  of 
mean-square-error) . 

To  go  beyond  Haar  wavelets,  new  spline  wavelet  refinement  operators  are  de¬ 
veloped.  These  operators  preserve  the  monotonicity  of  the  refinement. 
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2.  S+WAVELETS:  Algorithms  and  Technical  Details.  Andrew  Bruce, 
Hong- Ye  Gao,  and  David  Ragozin. 

A  complete  description  is  given  for  the  algorithms  in  S-|- WAVELETS  software 
toolkit  for  wavelet  analysis.  These  algorithms  include  wavelet  transforms, 
wavelet  packet  transforms,  cosine  packet  transforms,  and  non-decimated  wavelet 
transforms.  Implementations  for  the  transforms  and  their  inverses  are  given 
for  a  variety  of  boundary  treatment  rules,  including  periodic,  reflection,  in¬ 
terval  wavelets  (Cohen  et  al.  [CDV93]),  and  zero/polynomial  extension.  In 
addition,  modifications  to  the  standard  algorithms  are  given  to  handle  signals 
or  images  with  dimensions  not  divisible  by  a  power  of  two. 

3.5  Talks  and  Presentations 

1.  SPIE  Meetings.  Denoising  and  Robust  Non-Linear  Wavelet  Analysis,  April, 
1994,  Orlando,  FL. 

2.  CompStat  Meeting,  Smoothing  and  Robust  Wavelet  Analysis,  August,  1994, 
Vienna,  Austria. 

3.  IMS  Talks.  An  Object-Oriented  Toolkit  for  Wavelet  Analysis,  April,  1994, 
Cleveland,  OH. 

3.6  Papers  In  Preperation 

1.  Wavelet  Packet  Transforms  for  Finite  Signals.  Andrew  Bruce,  Hong- Ye 
Gao,  and  David  Ragozin. 

As  originally  developed  by  [Dau88]  the  discrete  wavelet  transform  is  defined 
on  the  entire  real  line.  To  apply  the  wavelet  and  wavelet  packet  transforms  to 
a  finite  signal,  several  approaches  have  been  pursued: 

(a)  An  ad  hoc  rule  is  used  for  applying  the  wavelet  filters  at  the  boundaries  of 
the  signal,  typically  by  recursively  extending  the  signal  at  each  iteration 
in  the  wavelet  decomposition. 

(b)  The  signal  is  extended  infinitely  according  to  some  boundary  rule  (e.g., 
periodic  extension)  and  sufficient  coefficients  are  retained  to  reproduce 
the  transform  (see,  for  example,  [MT92,  BF92]). 

(c)  A  new  wavelet  transform  is  formally  defined  on  compact  sets  (e.g.,  [CDV93] 
define  an  orthogonal  wavelet  basis  on  compact  intervals). 

A  unified  framework  is  developed  for  these  approaches.  The  methods  are  com¬ 
pared,  both  theoretically  and  empirically  with  emphasis  on  their  implications 
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for  data  analysis  and  statistical  estimation.  Practical  algorithms  are  given  for 
software  implementation. 

2.  A  Fast,  Robust,  Nonlinear  Wavelet  Transform.  Andrew  Bruce  and 
David  L.  Donoho. 

A  triadic  nonlinear  refinement  scheme  is  described  based  on  solving  the  fol¬ 
lowing  problem:  given  the  median  value  of  a  function  on  blocks  of  size  3~-', 
impute  medians  to  all  finer-scale  blocks  of  size  3“-'“'  by  fitting  local  polynomi¬ 
als  to  the  neighboring  block  medians  and  calculating  the  block  medians  of  the 
polynomials  one  scale  finer.  This  nonlinear  refinement  scheme  is  practical  for 
linear  and  quadratic  polynomials.  It  is  a  nonlinear  cousin  of  Deslauriers-Dubuc 
interpolation  and  of  Average-Interpolation.  It  is  reproduces  polynomials  up 
to  the  degree  of  the  fit. 

This  refinement  scheme  leads  to  nonlinear  multiresolution  analyses.  The  wavelet 
coefficient  functionals  are  possibly  nonlinear  combinations  of  a  finitely  many 
block  medians  at  the  same  scale;  a  fast  nonlinear  algorithm  to  compute  all 
block  medians  in  0{nlog{n))  time  can  be  made  by  a  simple  merging  of  sorted 
lists.  Using  this,  one  gets  a  fast  nonlinear  wavelet  transform. 

An  advantage  of  the  transform  is  its  robustness  against  outliers  and  impulsive 
noise,  while  preserving  high  degree  approximation  quality. 

3.  Spectral  Density  Estimation  by  Wavelet  Shrinkage.  Hong- Ye  Gao. 

We  study  the  problem  of  estimating  the  spectral  density  of  a  stationary  Gaus¬ 
sian  time  series.  We  use  an  orthogonal  wavelet  system  whose  members  are 
periodic  functions  and  have  a  finite  number  of  non-zero  Fourier  coefficients  - 
periodized  Meyer  wavelets.  We  apply  shrinkage  rules  to  the  empirical  wavelet 
coefficients.  We  show  that  estimates  based  on  thresholds  =  Xj  log  n  for 
certain  Aj,  with  n  the  sample  size,  have  near-optimal  convergence  rates, 
over  any  Besov  class  in  a  wide  range.  In  some  cases,  which  includes  the  Bump 
Algebra,  wavelet  shrinkage  procedures  significantly  outperform  classical  linear 
procedures,  such  as  window  methods  and  AR  approximation  methods. 

4.  8+ WAVELETS:  Applications  I  -  Wavelet  De-Noising.  Andrew  Bruce 

,  and  Hong-Ye  Gao. 

This  is  the  first  in  a  series  of  papers  illustrating  applications  of  wavelets  using 
S-f  Wavelets.  This  paper  focuses  on  the  removal  of  noise  in  a  variety  of 
applications.  Some  of  the  problems  examined  include  noise  removal  from  a 
noisy  NMR  signal,  signal  extraction  from  radar  glint  noise,  variance  estimate 
for  wavelet  de-noising,  image  enhancement,  probability  and  spectral  density 
estimation,  enhancement  of  noisy  speech  data,  and  isotonic  regression. 
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5.  S+WAVELETS:  Applications  II  -  Data  Compression.  Andrew  Bruce 
and  Hong- Ye  Gao. 

This  is  the  second  in  a  series  of  papers  illustrating  applications  of  wavelets 
using  S-h Wavelets.  Several  applications  of  wavelets  to  problems  in  data 
compression  are  presented.  Some  of  the  problems  include  compression  of  seis¬ 
mic  signals,  acoustic  signals,  digital  MR  images,  and  digital  fingerprint  images. 

6.  S+WAVELETS:  Applications  III  -  Signal  Processing.  Andrew  Bruce 
and  Hong- Ye  Gao. 

This  is  the  third  in  a  series  of  papers  illustrating  applications  of  wavelets 
using  S+ Wavelets.  A  variety  of  signal  processing  problems  is  examined, 
including  time-frequency  analysis  of  acoustic  signals,  singularity  detection  and 
estimation,  and  fast  algorithms  for  signal  classification. 

4  FUTURE  RESEARCH  PLANS 

Our  future  research  plans  will  continue  the  six  research  areas  as  discussed  in  sec¬ 
tion  2.  We  will  also  complete  the  papers  in  preparation  discussed  in  section  3.  In 

addition,  as  time  permits,  we  will  pursue  some  new  ideas  in  these  areas,  including 

Fast  Li/M-estimators  for  wavelets:  In  the  first  year  of  our  research,  we  found  that 
using  more  robust  methods  for  fitting  wavelets  produced  valuable  decomposi¬ 
tions  for  many  types  of  signals  and  images.  However,  these  methods  were  of 
limited  value  since  they  required  extensive  computational  effort.  Using  well 
known  stochastic  approximation  formula,  we  propose  to  develop  approximate 
methods  for  fast  computation  of  Li/M-estimators  for  wavelet  decompositions. 

Robust  time  series  filters.  The  ARMA  interpretation  of  the  wavelet  filters  will  be 
exploited  to  yield  a  robust  wavelet  decomposition  using  the  approximate  con¬ 
dition  mean  (ACM)  smoothers  developed  by  C.  J.  Masreliez  [Mas75]  and  R. 

D.  Martin  [Mar79,  MF93]. 

Robust  wavelets  in  higher  dimensions.  The  robust  algorithms  developed  for  one-dimensional 
signals  will  be  extended  to  higher-dimensional  data,  following  the  work  of 
[DJL92,  Ric93].  One  and  two  dimensional  non-decimated  algorithms  for  ro¬ 
bust  wavelet  analysis  will  also  be  explored. 

Wavelet  coefficients  for  fractional  Brownian  motion.  There  has  been  considerable 
activity  examining  the  statistical  properties  and  utility  of  wavelet  coefficients 
for  fractional  Brownian  motion.  In  particular,  there  is  work  on  texture  analysis 
[TD90],  estimation  and  analysis  of  fractional  Brownian  motion  [DT91,  Fla92], 
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extracting  fractal  signals  from  noisy  measurements  [W092],  correlation  struc¬ 
ture  of  wavelet  coefficients  for  fractional  Brownian  motion  [TK92,  DM94], 
Karhunen-Loeve  like  expansions  for  1//  processes  [Wor92],  and  analysis  of 
long-memory  processes  [PG94].  These  results  give  some  promise  for  further 
statistical  results. 

Simulation,  bootstrapping,  and  modeling  of  non- Gaussian  processes.  Simulation  of  non- 
Gaussian  processes  has  many  practical  uses  in  Navy  applications.  The  use 
of  wavelets  gives  hope  of  providing  a  convenient  method  for  simulating  non- 
Gaussian  random  processes.  One  would  compute  the  wavelet  coefficients  for 
an  observed  sample  of  a  non-Gaussian  process,  and  then  generate  simulation 
sample  paths  by  appropriate  random  sampling  of  the  wavelet  coefficients  and 
subsequent  expansion  in  terms  of  wavelets.  Various  methods  of  random  sam¬ 
pling  need  to  be  studied.  Then  one  needs  to  establish  that  the  sample  paths 
generated  converge  in  a  suitable  probabilistic  sense  to  the  probability  mea¬ 
sure  for  the  observed  sample  path  used  to  calculate  the  wavelet  coefficients  as 
the  sample  size  tends  to  infinity.  Work  in  this  area  could  lead  to  an  effective 
bootstrap  method  for  carrying  out  inference  for  non-Gaussian  processes.  This 
will  be  applied  to  improve  the  variance  estimation  methods  for  WaveShrink 
estimators. 
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Figure  1:  WaveShrink  with  non-decimated  wavelets.  Top:  the  noisy  jumpsine  signal. 
Bottom  (solid  line):  smooth  produced  by  wavelet  shrinkage  using  the  orthogonal 
DWT.  Bottom  (dashed  line):  smooth  produced  by  wavelet  shrinkage  using  the 
noil-decimated  DWT.  The  haar  wavelet  is  used  in  this  example. 
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Figure  2:  Approximate  confidence  band  for  WaveShrink.  Top:  the  annual  Nile  River 
flow  at  Aswan  from  1875  to  1970.  Middle:  WaveShrink  using  the  Haar  wavelet  with 
an  approximate  95%  confidence  band.  Bottom:  WaveShrink  using  the  bsl .  3  wavelet 
with  2-11  approximate  95%  confidence  band.  The  confidence  interval  clearly  shows  a 
significant  drop  around  the  time  of  the  building  of  the  Aswan  dam  in  1899-1902. 
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Figure  3:  First  row:  the  doppler  signal  (left)  and  the  noisy  doppler  signal  (right).  Second 
row:  the  WaveShrink  estimates  using  the  “universal”  threshold  (left)  and  the  “minimax” 
threshold  (right).  Third  row:  the  mean-square  error  (MSE)  of  the  WaveShrink  estimate 
of  every  other  sample  value  based  on  the  remaining  half  of  the  data  for  a  range  of  thresh¬ 
olds  from  1  to  3  (left)  and  the  cross-validation  WaveShrink  estimate  which  achieves  the 
minimum  MSE  (right).  Fourth  row:  MSE  between  the  WaveShrink  estimate  and  the  true 
doppler  signal  using  thresholds  ranging  from  1  to  3  (left)  and  the  WaveShrink  estimate  us¬ 
ing  the  threshold  which  achieves  the  minimum  (right).  Cross  validation  consistenly  leads 
to  lower  MSE  than  the  “universal”  and  “minimax”  thresholds  for  this  example. 
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Figure  4:  Comparison  of  wavelet  estimate  and  Grenander  estimate  for  monotone 
signal.  Top  left:  a  decreasing  curve  with  two  discontinuities,  Top  right:  the  same 
curve  plus  Gaussian  white  noise,  Bottom  left:  the  Grenander  estimate  with  mean- 
square  error.  Bottom  right:  the  wavelet  estimate  with  mean-square  error.  The 
wa■^elet  estimate  has  lower  mean  square  error  and  does  a  better  job  of  preserving 
the  jumps. 
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Figure  5:  A  simulated  AR  series  which  has  sharp  peaks  in  its  spectral  density 
function. 
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Figure  6:  The  true  spectrum  for  the  simulated  AR  series  of  figure  5  (top),  the  raw 
log-periodogram  (second),  the  smoothed  log-periodogram  using  a  triangular  spectral 
wiijdow  (third),  and  the  WaveShrink  estimate  (bottom). 
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Figure  7:  The  bumps  signal  (top),  the  bumps  signal  plus  Gaussian  white  noise 
(middle),  and  an  estimated  of  the  bumps  signal  obtained  from  WaveShrink  (bottom). 
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Figure  8:  Visualizing  the  WaveShrink  fit  for  the  bumps  signal.  Top  left:  the  de¬ 
composition  of  the  data  into  signal  plus  noise.  Top  right:  Boxplots  of  the  DWT 
coefficients  for  the  original  data  with  the  WaveShrink  thresholds  superimposed  as 
horizontal  lines.  Any  wavelet  coefficient  lying  inside  the  lines  is  set  to  zero.  Bottom 
left:  The  DWT  of  the  signal.  Bottom  right:  A  barplot  showing  the  decomposition 
of  the  energy  of  the  data  into  the  energy  attributable  to  signal  and  residual  energy. 
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Figure  9:  Four  views  of  the  residuals  from  the  WaveShrink  estimate  of  the  bumps 
signal.  Top  left:  the  DWT  of  the  residual  component.  Top  right:  the  autocorrelation 
function  (acf)  of  the  residual  component.  Bottom  left:  the  quantile-quantile  plot  of 
the^  residual  component  versus  the  quantiles  of  a  standard  normal.  Bottom  right:  a 
histogram  and  density  estimate  of  the  residuals. 
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Figure  10:  The  robust  smoother  algorithm  produces  a  pyramid  decomposition  with 
an  extra  component:  the  robust  residual  R{1).  For  each  multiresolution  level,  the 
low-pa^s  coefficients  S{t}  are  first  cleaned  using  a  robust  smoother  cleaner,  denoted 
by  sc  in  the  figure.  The  residuals  are  saved  in  the  R{1).  The  usual  wavelet  filters 
are  then  applied  to  the  cleaned  S{i)  to  obtain  S{1  +  1)  and  D{i  +  1). 


Figure  11:  (a)  Radar  glint  noise  data  in  degrees,  and  (b)  denoising  by  linear  shrinkage 
of  wavelets  (dashed  line)  compared  with  denoising  by  WaveShrink  combined  with 
robust  smoother-cleaner  wavelets  (solid  line). 
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Figure  12:  Three  acoustic  signals:  whale  clicks  (top),  shrimp  clicks  (middle),  and 
background  noise  (bottom). 
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Figure  13:  A  classification  tree  for  the  acoustic  signals  derived  from  a  wavelet  packet 
energy  map. 
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Figure  14:  A  wavelet  analysis  starts  with  the  transformation  of  a  signal  or  image 
object  into  a  wavelet  object.  A  wavelet  representation  is  not  a  single  transforma¬ 
tion,  but  involves  a  wide  range  of  objects.  The  toolkit  provides  a  rich  infrastructure 
for  interacting  with  the  data  and  wavelet  objects.  The  user  can  visualize,  manip¬ 
ulate  and  synthesize  wavelet  objects  with  built-in  functions  and  operators.  Be¬ 
cause  S-|- Wavelets  is  embedded  in  S-PLUS,  users  can  also  apply  a  wide  collection 
of  statistical  analysis  tools,  such  as  tree  network  classification,  linear  discriminant 
analysis,  or  non-linear  additive  models. 
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Figure  15:  An  overview  of  the  object-oriented  design  of  S-f- WAVELETS.  The  hori¬ 
zontal  boxes  group  the  objects  by  type:  1-D  transforms,  2-D  transforms,  etc..  The 
objects  are  also  categorized  as  primitive  objects  (middle  vertical  panel),  objects 
specialized  for  wavelet  and  wavelet  packet  analysis  (left  vertical  panel),  and  objects 
specialized  for  cosine  packet  analysis  (right  vertical  panel).  The  lines  indicate  the 
inheritance  hierarchy. 
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